{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First stab at visualizing ATL06\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import h5py\n",
    "import glob\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "# from bokeh.plotting import figure, output_file, show\n",
    "# from bokeh.models import ColumnDataSource, LinearColorMapper\n",
    "import numpy as np\n",
    "from ipyleaflet import Map, basemaps, basemap_to_tiles, GeoData, LayersControl\n",
    "import geopandas\n",
    "import json\n",
    "from readers.read_HDF5_ATL03 import read_HDF5_ATL03\n",
    "from readers.get_ATL03_x_atc import get_ATL03_x_atc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import getpass\n",
    "\n",
    "def list_files_local(path):\n",
    "    \"\"\" Get file list form local folder. \"\"\"\n",
    "    from glob import glob\n",
    "    return glob(path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<HDF5 file \"processed_ATL06_20190114020111_02540205_001_01.h5\" (mode r)>\n"
     ]
    }
   ],
   "source": [
    "prodname = 'ATL06'\n",
    "# prodname = 'ATL03'\n",
    "site = 'seeps-small-helheim'\n",
    "path2files = str('/home/jovyan/' + prodname + '/' + site + '/')\n",
    "# print(path2files)\n",
    "hdf5_flist = glob.glob(str(path2files +'*.h5'))\n",
    "# print(hdf5_flist)\n",
    "f = h5py.File(hdf5_flist[0], 'r')\n",
    "print(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyproj\n",
    "from astropy.time import Time\n",
    "\n",
    "def gps2dyr(time):\n",
    "    \"\"\" Converte GPS time to decimal years. \"\"\"\n",
    "    return Time(time, format='gps').decimalyear\n",
    "\n",
    "def track_type(time, lat, tmax=1):\n",
    "    \"\"\"\n",
    "    Separate tracks into ascending and descending.\n",
    "    \n",
    "    Defines tracks as segments with time breaks > tmax,\n",
    "    and tests whether lat increases or decreases w/time.\n",
    "    \"\"\"\n",
    "    tracks = np.zeros(lat.shape)  # generate track segment\n",
    "    tracks[0:np.argmax(np.abs(lat))] = 1  # set values for segment\n",
    "    i_asc = np.zeros(tracks.shape, dtype=bool)  # output index array\n",
    "\n",
    "    # Loop trough individual secments\n",
    "    for track in np.unique(tracks):\n",
    "    \n",
    "        i_track, = np.where(track == tracks)  # get all pts from seg\n",
    "    \n",
    "        if len(i_track) < 2: continue\n",
    "    \n",
    "        # Test if lat increases (asc) or decreases (des) w/time\n",
    "        i_min = time[i_track].argmin()\n",
    "        i_max = time[i_track].argmax()\n",
    "        lat_diff = lat[i_track][i_max] - lat[i_track][i_min]\n",
    "    \n",
    "        # Determine track type\n",
    "        if lat_diff > 0:  i_asc[i_track] = True\n",
    "    \n",
    "    return i_asc, np.invert(i_asc)  # index vectors\n",
    "\n",
    "\n",
    "def transform_coord(proj1, proj2, x, y):\n",
    "    \"\"\"\n",
    "    Transform coordinates from proj1 to proj2 (EPSG num).\n",
    "\n",
    "    Example EPSG projs:\n",
    "        Geodetic (lon/lat): 4326\n",
    "        Polar Stereo AnIS (x/y): 3031\n",
    "        Polar Stereo GrIS (x/y): 3413\n",
    "    \"\"\"\n",
    "    # Set full EPSG projection strings\n",
    "    proj1 = pyproj.Proj(\"+init=EPSG:\"+str(proj1))\n",
    "    proj2 = pyproj.Proj(\"+init=EPSG:\"+str(proj2))\n",
    "    return pyproj.transform(proj1, proj2, x, y)  # convert\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib widget\n",
    "\n",
    "def read_h5(fname, vnames=[]):\n",
    "    \"\"\" Simple HDF5 reader. \"\"\"\n",
    "    with h5py.File(fname, 'r') as f:\n",
    "        return [f[v][:] for v in vnames]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190125125128_04290203_001_01.h5\n"
     ]
    }
   ],
   "source": [
    "# files = list_files_local(path2files)\n",
    "# print(files)\n",
    "fname = hdf5_flist[1]\n",
    "print(fname)\n",
    "\n",
    "#!h5ls -r /home/jovyan/ATL06/seeps/processed_ATL06_20181030170315_04900103_001_01.h5\n",
    "# !h5ls -r /home/jovyan/ATL03/seeps/processed_ATL03_20190202123455_05510203_001_01.h5\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "\n",
    "def read_atl06(fname, bbox=None):\n",
    "    \"\"\" \n",
    "    Read 1 ATL06 file and output 6 reduced files. \n",
    "    \n",
    "    Extract variables of interest and separate the ATL06 file \n",
    "    into each beam (ground track) and ascending/descending orbits.\n",
    "    \"\"\"\n",
    "\n",
    "    # Each beam is a group\n",
    "    group = ['/gt1l', '/gt1r', '/gt2l', '/gt2r', '/gt3l', '/gt3r']\n",
    "    # Loop trough beams\n",
    "    for k,g in enumerate(group):\n",
    "\n",
    "        #-----------------------------------#\n",
    "        # 1) Read in data for a single beam #\n",
    "        #-----------------------------------#\n",
    "    \n",
    "        # Load variables into memory (more can be added!)\n",
    "        with h5py.File(fname, 'r') as fi:\n",
    "            #print(fi[g].keys())\n",
    "            if g+'/land_ice_segments' in fi.keys():\n",
    "                #print('Found some data on beam' + g)\n",
    "                #print(fi[g+'/land_ice_segments'].keys())                \n",
    "\n",
    "                lat = fi[g+'/land_ice_segments/latitude'][:]\n",
    "                lon = fi[g+'/land_ice_segments/longitude'][:]\n",
    "                h_li = fi[g+'/land_ice_segments/h_li'][:]\n",
    "                s_li = fi[g+'/land_ice_segments/h_li_sigma'][:]\n",
    "                t_dt = fi[g+'/land_ice_segments/delta_time'][:]\n",
    "                q_flag = fi[g+'/land_ice_segments/atl06_quality_summary'][:]\n",
    "                t_ref = fi['/ancillary_data/atlas_sdp_gps_epoch'][:]\n",
    "                rgt = fi['/orbit_info/rgt'][:] * np.ones(len(lat))\n",
    "                orb = np.full_like(h_li, k)\n",
    "\n",
    "                #---------------------------------------------#\n",
    "                # 2) Filter data according region and quality #\n",
    "                #---------------------------------------------#\n",
    "\n",
    "                # Select a region of interest\n",
    "                if bbox:\n",
    "                    lonmin, lonmax, latmin, latmax = bbox\n",
    "                    bbox_mask = (lon >= lonmin) & (lon <= lonmax) & \\\n",
    "                                (lat >= latmin) & (lat <= latmax)\n",
    "                else:\n",
    "                    bbox_mask = np.ones_like(lat, dtype=bool)  # get all\n",
    "\n",
    "                # Only keep good data, and data inside bbox\n",
    "                mask = (q_flag == 0) & (np.abs(h_li) < 10e3) & (bbox_mask == 1)\n",
    "\n",
    "                # Update variables\n",
    "                lat, lon, h_li, s_li, t_dt, rgt, orb = \\\n",
    "                    lat[mask], lon[mask], h_li[mask], s_li[mask], t_dt[mask],\\\n",
    "                        rgt[mask], orb[mask]\n",
    "\n",
    "                # Test for no data\n",
    "                if len(h_li) == 0: continue\n",
    "\n",
    "                #-------------------------------------#\n",
    "                # 3) Convert time and separate tracks #\n",
    "                #-------------------------------------#\n",
    "\n",
    "                # Time in GPS seconds (secs sinde 1980...)\n",
    "                t_gps = t_ref + t_dt\n",
    "\n",
    "                # Time in decimal years\n",
    "                t_year = gps2dyr(t_gps)\n",
    "\n",
    "                # Determine orbit type\n",
    "                i_asc, i_des = track_type(t_year, lat)\n",
    "\n",
    "                #-----------------------#\n",
    "                # 4) Save selected data #\n",
    "                #-----------------------#\n",
    "\n",
    "                # Define output file name\n",
    "                ofile = fname.replace('.h5', '_'+g[1:]+'.h5')\n",
    "\n",
    "                # Save variables\n",
    "                with h5py.File(ofile, 'w') as f:\n",
    "                    f['orbit'] = orb\n",
    "                    f['lon'] = lon\n",
    "                    f['lat'] = lat\n",
    "                    f['h_elv'] = h_li\n",
    "                    f['t_year'] = t_year\n",
    "                    f['t_sec'] = t_gps\n",
    "                    f['s_elv'] = s_li\n",
    "                    f['q_flg'] = q_flag\n",
    "                    f['rgt'] = rgt\n",
    "                    f['trk_type'] = i_asc\n",
    "\n",
    "                    print('out ->', ofile)\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "\n",
    "def read_atl03(fname, bbox=None):\n",
    "    print(fname)\n",
    "    \"\"\" \n",
    "    Read 1 ATL06 file and output 6 reduced files. \n",
    "    \n",
    "    Extract variables of interest and separate the ATL06 file \n",
    "    into each beam (ground track) and ascending/descending orbits.\n",
    "    \"\"\"\n",
    "    f = h5py.File(fname,'r')\n",
    "    print(f.keys())\n",
    "    # Each beam is a group\n",
    "    group = ['/gt1l', '/gt1r', '/gt2l', '/gt2r', '/gt3l', '/gt3r']\n",
    "    # Loop trough beams\n",
    "    for k,g in enumerate(group):\n",
    "    \n",
    "        #-----------------------------------#\n",
    "        # 1) Read in data for a single beam #\n",
    "        #-----------------------------------#\n",
    "    \n",
    "        # Load variables into memory (more can be added!)\n",
    "        with h5py.File(fname, 'r') as fi:\n",
    "            print(fi[g].keys())\n",
    "            if g+'/heights' in fi.keys():\n",
    "                \n",
    "                print('Found some data on beam' + g)\n",
    "                print(fi[g+'/heights'].keys()) \n",
    "                \n",
    "                lat = fi[g+'/heights/lat_ph'][:]\n",
    "                lon = fi[g+'/heights/lon_ph'][:]\n",
    "                h = fi[g+'/heights/h_ph'][:]\n",
    "                conf = fi[g+'/heights/signal_conf_ph']\n",
    "\n",
    "                land_ice_class = conf[:,3]\n",
    "          \n",
    "        #-----------------------------------#\n",
    "        # 3) Filter data #\n",
    "        #-----------------------------------#\n",
    "                mask = (land_ice_class == 4) & (np.abs(h) < 10e3)\n",
    "                lat,lon,h = lat[mask],lon[mask],h[mask]\n",
    "        \n",
    "        #-----------------------#\n",
    "        # 4) Save selected data #\n",
    "        #-----------------------#\n",
    "        \n",
    "        # Define output file name\n",
    "                ofile = fname.replace('.h5', '_'+g[1:]+'.h5')\n",
    "                \n",
    "        # Save variables\n",
    "                with h5py.File(ofile, 'w') as f:\n",
    "                    f['lon'] = lon\n",
    "                    f['lat'] = lat\n",
    "                    f['h_elv'] = h\n",
    "            \n",
    "                    print('out ->', ofile)\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file number 0 is /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01_gt1l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01_gt1r.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01_gt2l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01_gt2r.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01_gt3l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01_gt3r.h5\n",
      "file number 1 is /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190125125128_04290203_001_01.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190125125128_04290203_001_01_gt1l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190125125128_04290203_001_01_gt1r.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190125125128_04290203_001_01_gt2l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190125125128_04290203_001_01_gt2r.h5\n",
      "file number 2 is /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181026171133_04290103_001_01.h5\n",
      "file number 3 is /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190212003725_06960205_001_01.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190212003725_06960205_001_01_gt2l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190212003725_06960205_001_01_gt2r.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190212003725_06960205_001_01_gt3l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190212003725_06960205_001_01_gt3r.h5\n",
      "file number 4 is /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01_gt1l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01_gt1r.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01_gt2l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01_gt2r.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01_gt3l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01_gt3r.h5\n",
      "file number 5 is /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181128153912_09320103_001_01.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181128153912_09320103_001_01_gt2l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181128153912_09320103_001_01_gt2r.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181128153912_09320103_001_01_gt3l.h5\n",
      "out -> /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181128153912_09320103_001_01_gt3r.h5\n"
     ]
    }
   ],
   "source": [
    "# print(hdf5_flist)\n",
    "for ii in range(0,len(hdf5_flist)):\n",
    "    print('file number ' + str(ii) + ' is ' + hdf5_flist[ii])\n",
    "    read_atl06(hdf5_flist[ii])\n",
    "    #read_atl03(hdf5_flist[ii])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/land_ice_segments       Group\n",
      "/land_ice_segments/atl06_quality_summary Dataset {792/Inf}\n",
      "/land_ice_segments/delta_time Dataset {792/Inf}\n",
      "/land_ice_segments/h_li  Dataset {792/Inf}\n",
      "/land_ice_segments/h_li_sigma Dataset {792/Inf}\n",
      "/land_ice_segments/latitude Dataset {792/Inf}\n",
      "/land_ice_segments/longitude Dataset {792/Inf}\n",
      "/land_ice_segments/segment_id Dataset {792/Inf}\n",
      "/land_ice_segments/sigma_geo_h Dataset {792/Inf}\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hdf_files_by_laser = list_files_local(path2files + '*gt*.h5')\n",
    "#!h5ls -r /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190125125128_04290203_001_01_gt2r.h5\n",
    "# !h5ls -r /home/jovyan/ATL03/seeps-small-helheim/processed_ATL03_20190125125128_04290203_001_01_gt2r.h5\n",
    "\n",
    "\n",
    "#!h5ls -r /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01.h5/gt2l/\n",
    "# !h5ls -r /home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20190114020111_02540205_001_01_gt2l.h5\n",
    "\n",
    "\n",
    "%matplotlib widget\n",
    "plt.figure()\n",
    "#print(hdf_files_by_laser)\n",
    "for ii in range(0, len(hdf_files_by_laser)):\n",
    "    with h5py.File(hdf_files_by_laser[ii]) as f:\n",
    "        #print(f)\n",
    "        elev_atl06 = f['h_elv'][:]\n",
    "        sigma_atl06 = f['s_elv'][:]\n",
    "        lat_atl06 = f['/lat'][:]\n",
    "        lon_atl06 = f['/lon'][:]\n",
    "        #background_atl06 = f['/gt2l/land_ice_segments/geophysical/bckgrd'][:]\n",
    "        #ebackground_atl06 = f['/gt2l/land_ice_segments/geophysical/e_bckgrd'][:]\n",
    "        f.close()\n",
    "#         plt.scatter(lon_atl06, lat_atl06, c=sigma_atl06)\n",
    "        plt.scatter(lon_atl06, lat_atl06, c=elev_atl06)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "792\n",
      "792\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe63bd67978>]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Grabbing the data, sticking to the center strong laser\n",
    "\n",
    "with h5py.File('/home/jovyan/ATL06/seeps-small-helheim/processed_ATL06_20181227141518_13740103_001_01.h5') as f:\n",
    "    elev_land_ice = f['/gt2l/land_ice_segments/h_li'][:]\n",
    "    sigma_land_ice = f['gt2l/land_ice_segments/h_li_sigma'][:]\n",
    "    qual = f['gt2l/land_ice_segments/atl06_quality_summary'][:]\n",
    "#     xatc_land_ice = f['/gt2l/land_ice_segments/ground_track/x_atc'][:]\n",
    "\n",
    "print(len(sigma_land_ice))\n",
    "print(len(elev_land_ice))\n",
    "elev_land_ice[elev_land_ice > 10000] = np.nan\n",
    "\n",
    "%matplotlib widget\n",
    "plt.figure()\n",
    "plt.plot(elev_land_ice)\n",
    "plt.plot(elev_land_ice+sigma_land_ice)\n",
    "plt.plot(elev_land_ice-sigma_land_ice)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6c2008b677e6447e874b9e421a8eb95d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe625d17128>]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "plt.figure()\n",
    "plt.plot(sigma_land_ice)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
